### Import statewide pedestrian crash data 2017-2019 data and preview
pedcrashes<-read_csv("../Data/pedcrashes20172019.csv")
## Rows: 6079 Columns: 59
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (48): crash_numb, city_town_name, crash_date, crash_severity_descr, cra...
## dbl  (10): year, numb_vehc, district_num, numb_fatal_injr, numb_nonfatal_inj...
## time  (1): crash_time_2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
paged_table((pedcrashes))
## Warning: One or more parsing issues, see `problems()` for details
### Select fields of interest and preview
pedcrashes<-pedcrashes%>%
  select(crash_numb, max_injr_svrty_cl, drvr_cntrb_circ_cl, road_surf_cond_descr, road_cntrb_descr, traf_cntrl_devc_type_descr, ambnt_light_descr)
paged_table(head(pedcrashes))
###Summarize fields of interest
summary(pedcrashes)
##   crash_numb        max_injr_svrty_cl  drvr_cntrb_circ_cl road_surf_cond_descr
##  Length:6079        Length:6079        Length:6079        Length:6079         
##  Class :character   Class :character   Class :character   Class :character    
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character    
##  road_cntrb_descr   traf_cntrl_devc_type_descr ambnt_light_descr 
##  Length:6079        Length:6079                Length:6079       
##  Class :character   Class :character           Class :character  
##  Mode  :character   Mode  :character           Mode  :character
### Identify unique values of Max Injury Severity Reported
table(select(pedcrashes, max_injr_svrty_cl))
## 
##                      Fatal injury (K)                No Apparent Injury (O) 
##                                   117                                   236 
##                             No injury     Non-fatal injury - Incapacitating 
##                                   867                                   539 
## Non-fatal injury - Non-incapacitating           Non-fatal injury - Possible 
##                                  1804                                  1522 
##                   Possible Injury (C)            Suspected Minor Injury (B) 
##                                   265                                   552 
##          Suspected Serious Injury (A) 
##                                   174
##Use case_when to create new field "injury_level" and mutate max_injr_svrty_cl into numeric scale using the KABCO injury classification scale
#### K=5, A=4, B=3, C=2, O=1
pedcrashes1<-pedcrashes%>%
  mutate(injury_level = case_when(
         max_injr_svrty_cl == "Fatal injury (K)" ~ 5,
         max_injr_svrty_cl == "Non-fatal injury - Incapacitating"~ 4,
         max_injr_svrty_cl == "Suspected Serious Injury (A)" ~ 4,
         max_injr_svrty_cl == "Non-fatal injury - Non-incapacitating" ~ 3,
         max_injr_svrty_cl == "Suspected Minor Injury (B)" ~ 3,
         max_injr_svrty_cl == "Non-fatal injury - Possible" ~ 2,
         max_injr_svrty_cl == "Possible Injury (C)" ~ 2,
         max_injr_svrty_cl == "No injury" ~ 1,
         max_injr_svrty_cl == "No Apparent Injury (O)" ~ 1,
         TRUE ~ -1))
  
pedcrashesgraph<-pedcrashes1%>%
  filter(injury_level != -1) %>%
  mutate(injury_label = case_when(
         injury_level == 5 ~ "Fatal injury (K)",
         injury_level == 4 ~ "Suspected Serious Injury (A)",
         injury_level == 3 ~ "Possible Injury (C)",
         injury_level == 2 ~ "Suspected Minor Injury (B)",
         injury_level == 1 ~ "No Apparent Injury (O)",
         TRUE ~ "NA"))
### Create bar chart for crash severity distribution

barchart<-pedcrashesgraph%>%
  ggplot(aes(injury_label)) +geom_bar(stat="count",width=.7, fill="lightblue") +
  theme(axis.text.x = element_text(angle = 10, size= 13, vjust= 1, hjust=1))+ ### adjust x-axis label orientation
  geom_text(stat='count', aes(label=..count..), size= 6, vjust=1.6 )+ ### add count labels
  ggtitle("Massachusetts 2017-2019 Pedestrian Crashes by Crash Severity")+
  xlab("Crash Severity")+
  ylab("Number of Crashes")+
  theme(plot.title = element_text(size=22, face="bold"),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))

barchart

### Export bar chart

ggsave("pedcrashbarchart.png", plot=barchart,width = 330, 
       height = 200, 
       units = "mm")
##Use case_when to create new field "is_ka" and mutate max_injr_svrty_cl into dummy for 1= KA, 0 = not a KA crash

pedcrashes2<-pedcrashes1%>%
  mutate(is_ka = case_when(
         max_injr_svrty_cl == "Fatal injury (K)" ~ 1,
         max_injr_svrty_cl == "Non-fatal injury - Incapacitating"~ 1,
         max_injr_svrty_cl == "Suspected Serious Injury (A)" ~ 1,
         max_injr_svrty_cl == "Non-fatal injury - Non-incapacitating" ~ 0,
         max_injr_svrty_cl == "Suspected Minor Injury (B)" ~ 0,
         max_injr_svrty_cl == "Non-fatal injury - Possible" ~ 0,
         max_injr_svrty_cl == "Possible Injury (C)" ~ 0,
         max_injr_svrty_cl == "No injury" ~ 0,
         max_injr_svrty_cl == "No Apparent Injury (O)" ~ 0,
         TRUE ~ -1))

table(select(pedcrashes2, is_ka))
## 
##   -1    0    1 
##    3 5246  830
### Identify unique values of Road Surface Condition description
table(select(pedcrashes, road_surf_cond_descr))
## 
##                          Dry                          Ice 
##                         4925                           29 
## Sand, mud, dirt, oil, gravel                        Slush 
##                           12                           11 
##                         Snow     Water (standing, moving) 
##                           53                            1 
##                          Wet 
##                         1045
### create dummy road_surface where 1 = road surface altered, 0 = dry road surface
pedcrashes3<-pedcrashes2%>%
  mutate(road_surface = case_when(
         road_surf_cond_descr == "Dry" ~ 0,
         road_surf_cond_descr == "Ice"~ 1,
         road_surf_cond_descr == "Sand, mud, dirt, oil, gravel"~ 1,
         road_surf_cond_descr == "Slush" ~ 1,
         road_surf_cond_descr == "Snow" ~ 1,
         road_surf_cond_descr == "Water (standing, moving)" ~ 1,
         road_surf_cond_descr == "Wet" ~ 1,
         road_cntrb_descr == "Road surface condition (wet, icy, snow, slush, etc.)" ~ 1,
         TRUE ~ -1))

table(select(pedcrashes3, road_surface))
## 
##   -1    0    1 
##    3 4925 1151
### Identify unique values of Road Contributing Circumstance Descriptions
table(select(pedcrashes, road_cntrb_descr))
## 
##                                                   Debris 
##                                                        3 
##                                                     None 
##                                                     5565 
##                                   Obstruction in roadway 
##                                                       17 
##                                                    Other 
##                                                       79 
##     Road surface condition (wet, icy, snow, slush, etc.) 
##                                                      202 
##                                        Rut, holes, bumps 
##                                                        6 
##                              Shoulders (none, low, soft) 
##                                                        4 
##                                 Toll/booth/plaza related 
##                                                        2 
##                               Traffic congestion related 
##                                                      148 
## Traffic control device inoperative, missing, or obscured 
##                                                       13 
##             Work zone (construction/maintenance/utility) 
##                                                       32 
##                            Worn, travel-polished surface 
##                                                        5
### create dummy road_circumstances where 1 = None , 0 = road contributing circumstance
pedcrashes4<-pedcrashes3%>%
  mutate(road_circumstances = case_when(
         road_cntrb_descr == "None" ~ 0,
         road_cntrb_descr == "Debris"~ 1,
         road_cntrb_descr == "Obstruction in roadway"~ 1,
         road_cntrb_descr == "Other" ~ 1,
         road_cntrb_descr == "Road surface condition (wet, icy, snow, slush, etc.)" ~ 0, ###coded 0 since value accounted for in road_surface
         road_cntrb_descr == "Rut, holes, bumps" ~ 1,
         road_cntrb_descr == "Shoulders (none, low, soft)" ~ 1,
         road_cntrb_descr == "Toll/booth/plaza related" ~ 1,
         road_cntrb_descr == "Traffic congestion related" ~ 1,
         road_cntrb_descr == "Traffic control device inoperative, missing, or obscured" ~ 1,
         road_cntrb_descr == "Work zone (construction/maintenance/utility)" ~ 1,
         road_cntrb_descr == "Worn, travel-polished surface" ~ 1,
         TRUE ~ -1))

table(select(pedcrashes4, road_circumstances))
## 
##   -1    0    1 
##    3 5767  309
### Identify unique values of Ambient Light Descriptions
table(select(pedcrashes,ambnt_light_descr))
## 
##          Dark - lighted roadway      Dark - roadway not lighted 
##                            1402                             270 
## Dark - unknown roadway lighting                            Dawn 
##                              40                              71 
##                        Daylight                            Dusk 
##                            4130                             163
### create dummy ambient_light where 1 = Daylight, 0 = Dark or potential lack of light
pedcrashes5<-pedcrashes4%>%
  mutate(ambient_light = case_when(
         ambnt_light_descr == "Dark - lighted roadway" ~ 0,
         ambnt_light_descr == "Dark - unknown roadway lighting"~ 0,
         ambnt_light_descr == "Dark - roadway not lighted" ~ 0,
         ambnt_light_descr == "Dawn" ~ 0,
         ambnt_light_descr == "Daylight" ~ 1,
         ambnt_light_descr == "Dusk" ~ 0,
         TRUE ~ -1))

table(select(pedcrashes5, ambient_light))
## 
##   -1    0    1 
##    3 1946 4130
### Identify unique values of Road Contributing Circumstance Descriptions
table(select(pedcrashes, traf_cntrl_devc_type_descr))
## 
## Flashing traffic control signal                     No controls 
##                              65                            3807 
##         Railway crossing device               School zone signs 
##                               2                              29 
##                      Stop signs          Traffic control signal 
##                             856                            1162 
##                   Warning signs                     Yield signs 
##                              79                              76
### create dummy traffic_control_device where 1 = no traffic control device present, 0 = no traffic controls present
pedcrashes6<-pedcrashes5%>%
  mutate(traffic_control_device = case_when(
         traf_cntrl_devc_type_descr == "No controls" ~ 0,
         traf_cntrl_devc_type_descr == "Yield signs"~ 1,
         traf_cntrl_devc_type_descr == "Traffic control signal"~ 1,
         traf_cntrl_devc_type_descr == "Flashing traffic control signal" ~ 1,
         traf_cntrl_devc_type_descr == "Warning signs" ~ 1,
         traf_cntrl_devc_type_descr == "Railway crossing device" ~ 1,
         traf_cntrl_devc_type_descr == "Yield signs" ~ 1,
         traf_cntrl_devc_type_descr == "School zone signs" ~ 1,
         traf_cntrl_devc_type_descr == "Stop signs" ~ 1,
         TRUE ~ -1))

table(select(pedcrashes6, traffic_control_device))
## 
##   -1    0    1 
##    3 3807 2269
### Identify unique values of Driver Contributing Circumstance Descriptions
## Many values, remove comment to run: table(select(pedcrashes, drvr_cntrb_circ_cl))
### create dummy traffic_control_device where 1 = at least 1 driver had improper driving, 0 = no improper driving for at least 1 driver, others may be unknown 
pedcrashes7<-pedcrashes6%>%
  mutate(improper_driving =  case_when(
    str_detect(drvr_cntrb_circ_cl, "speed limit") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "yield") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "too closely") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "improper turn") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "visibility") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Visibility") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "too fast") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Distracted") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Inattention") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "reckless") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "aggressive") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "steering") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Disregarded traffic signs") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Other improper action") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Swerving") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Glare") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Failure") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Illness") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Fatigue") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Emotional") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "History heart/epilepsy/fainting") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Wrong side or wrong way") ~ 1,
    str_detect(drvr_cntrb_circ_cl, "Physical impairment") ~ 1,
    TRUE ~ 0))

table(select(pedcrashes7, improper_driving))
## 
##    0    1 
## 3134 2945
### Check if all driver contributing circumstances are captured, excluding defective equipment
filtercheck<-pedcrashes7%>%
  filter(improper_driving == 0) %>%
  paged_table()
table(select(filtercheck, drvr_cntrb_circ_cl))
## 
##                                                                                     D1: (),(No improper driving) 
##                                                                                                                7 
##                                                                                        D1: (No improper driving) 
##                                                                                                             2937 
##                                                           D1: (No improper driving)  / D2: (No improper driving) 
##                                                                                                               28 
##                              D1: (No improper driving)  / D2: (No improper driving)  / D3: (No improper driving) 
##                                                                                                                2 
## D1: (No improper driving)  / D2: (No improper driving)  / D3: (No improper driving)  / D4: (No improper driving) 
##                                                                                                                1 
##                                     D1: (No improper driving)  / D2: (No improper driving),(No improper driving) 
##                                                                                                                1 
##                                                                       D1: (No improper driving)  / D2: (Unknown) 
##                                                                                                                3 
##                                                           D1: (No improper driving)  / D4: (No improper driving) 
##                                                                                                                1 
##                                                                  D1: (No improper driving),(No improper driving) 
##                                                                                                              120 
##                                                                              D1: (No improper driving),(Unknown) 
##                                                                                                                9 
##                                                                              D1: (Operating defective equipment) 
##                                                                                                                2 
##                                                                       D1: (Unknown)  / D2: (No improper driving) 
##                                                                                                                5 
##                                          D1: (Unknown)  / D2: (No improper driving)  / D3: (No improper driving) 
##                                                                                                                1 
##                                                                       D1: (Unknown)  / D3: (No improper driving) 
##                                                                                                                1 
##                                                      D1: (Unknown)  / D3: (No improper driving)  / D4: (Unknown) 
##                                                                                                                1 
##                                                                              D1: (Unknown),(No improper driving) 
##                                                                                                                5 
##                                       D1: (Unknown),(Unknown)  / D2: (No improper driving),(No improper driving) 
##                                                                                                                1 
##                                                                                        D2: (No improper driving) 
##                                                                                                                4 
##                                                           D2: (No improper driving)  / D3: (No improper driving) 
##                                                                                                                1 
##               D3: (No improper driving),(No improper driving)  / D4: (No improper driving),(No improper driving) 
##                                                                                                                1
##Filter to only necessary fields
###Filter out unknowns for all variables and view data with relevant fields
pedcrashesclean<-pedcrashes7%>%
  select(injury_level, is_ka, road_surface, road_circumstances, ambient_light, traffic_control_device, improper_driving)%>%
  filter(traffic_control_device!= -1, na.rm= TRUE)

summary(pedcrashesclean)
##   injury_level       is_ka         road_surface    road_circumstances
##  Min.   :1.000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   
##  1st Qu.:2.000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000   
##  Median :3.000   Median :0.0000   Median :0.0000   Median :0.00000   
##  Mean   :2.499   Mean   :0.1366   Mean   :0.1894   Mean   :0.05086   
##  3rd Qu.:3.000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.00000   
##  Max.   :5.000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000   
##  ambient_light    traffic_control_device improper_driving
##  Min.   :0.0000   Min.   :0.0000         Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000         1st Qu.:0.0000  
##  Median :1.0000   Median :0.0000         Median :0.0000  
##  Mean   :0.6797   Mean   :0.3734         Mean   :0.4847  
##  3rd Qu.:1.0000   3rd Qu.:1.0000         3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :1.0000         Max.   :1.0000
###Create Training data
set.seed(1)
trainingsetsize<-pedcrashesclean%>%
  nrow()%>%
  multiply_by(0.7)%>%
  round()
train_indices<-sample(1:nrow(pedcrashesclean), trainingsetsize)
train<-pedcrashesclean[train_indices,]
test<-pedcrashesclean[-train_indices,]
ols_levels<-lm(injury_level ~ .-is_ka, data = train)
ols_binary<-lm(is_ka ~ .-injury_level, data = train)
par(mfrow = c(2,3))
plot(ols_binary, which = 1:6)

par(mfrow = c(2,3))
plot(ols_levels, which = 1:6)

# backward elimination
step(object = ols_levels, # specify full model
     direction = "backward") # specify backward
## Start:  AIC=-247.45
## injury_level ~ (is_ka + road_surface + road_circumstances + ambient_light + 
##     traffic_control_device + improper_driving) - is_ka
## 
##                          Df Sum of Sq    RSS     AIC
## - road_surface            1    1.1283 4002.4 -248.25
## - road_circumstances      1    1.3730 4002.7 -247.99
## <none>                                4001.3 -247.45
## - traffic_control_device  1    9.2060 4010.5 -239.68
## - improper_driving        1   24.6577 4026.0 -223.32
## - ambient_light           1   25.5457 4026.9 -222.38
## 
## Step:  AIC=-248.25
## injury_level ~ road_circumstances + ambient_light + traffic_control_device + 
##     improper_driving
## 
##                          Df Sum of Sq    RSS     AIC
## - road_circumstances      1    1.3532 4003.8 -248.81
## <none>                                4002.4 -248.25
## - traffic_control_device  1    9.3477 4011.8 -240.33
## - ambient_light           1   24.4184 4026.9 -224.38
## - improper_driving        1   24.4951 4026.9 -224.30
## 
## Step:  AIC=-248.81
## injury_level ~ ambient_light + traffic_control_device + improper_driving
## 
##                          Df Sum of Sq    RSS     AIC
## <none>                                4003.8 -248.81
## - traffic_control_device  1    9.4998 4013.3 -240.74
## - improper_driving        1   24.6447 4028.4 -224.72
## - ambient_light           1   24.8112 4028.6 -224.54
## 
## Call:
## lm(formula = injury_level ~ ambient_light + traffic_control_device + 
##     improper_driving, data = train)
## 
## Coefficients:
##            (Intercept)           ambient_light  traffic_control_device  
##                2.58184                -0.16408                -0.09745  
##       improper_driving  
##                0.15285
###test backwards elimination model
br_fit<-lm(formula = injury_level ~ ambient_light + traffic_control_device + improper_driving, data = train)
summary(br_fit)
## 
## Call:
## lm(formula = injury_level ~ ambient_light + traffic_control_device + 
##     improper_driving, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7347 -0.5818  0.2653  0.5822  2.5822 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             2.58184    0.03118  82.797  < 2e-16 ***
## ambient_light          -0.16408    0.03198  -5.131 3.01e-07 ***
## traffic_control_device -0.09745    0.03069  -3.175  0.00151 ** 
## improper_driving        0.15285    0.02989   5.114 3.29e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9707 on 4249 degrees of freedom
## Multiple R-squared:  0.0136, Adjusted R-squared:  0.01291 
## F-statistic: 19.53 on 3 and 4249 DF,  p-value: 1.422e-12
###Create and print diagnostics plots for backwards regression model
png("diagnosticsplots.png", width=15, height=6, units='in', res=300)
par(mfrow = c(1:2))
plot(br_fit, which=1:2, ask=FALSE)
dev.off()
## quartz_off_screen 
##                 2
###Model Evaluation- run models using testing data, generate AIC, BIC, PRESS, MSE, RSME

predict(br_fit, test) %>% subtract(test$injury_level) %>%
  raise_to_power(2) %>% mean() -> MSE_br_fit
predict(ols_binary, test) %>% subtract(test$is_ka) %>%
  raise_to_power(2) %>% mean() -> MSE_ols_binary
predict(ols_levels, test) %>% subtract(test$injury_level) %>% 
  raise_to_power(2) %>% mean() -> MSE_ols_levels
Press1<-press(br_fit)
BIC1<-BIC(br_fit)
AIC1<-AIC(br_fit)
Press2<-press(ols_levels)
BIC2<-BIC(ols_levels)
AIC2<-AIC(ols_levels)
###visualize metrics in table
cat('               ', " MSE", '      RMSE', '     PRESS','   AIC', '     BIC', '\n',
'Backwards Elim:', MSE_br_fit, sqrt(MSE_br_fit), Press1, AIC1, BIC1, '\n',
'All Predictors:', MSE_ols_levels, sqrt(MSE_ols_levels), Press2, AIC2,BIC2, '\n')
##                  MSE       RMSE      PRESS    AIC      BIC 
##  Backwards Elim: 0.9702171 0.984996 4011.383 11822.68 11854.45 
##  All Predictors: 0.9714831 0.9856384 4012.704 11824.04 11868.53
stargazer(ols_levels, br_fit, title="Table 1. Regression Results",
dep.var.labels=c("Injury Severity Level"),
covariate.labels=c("Pavement Friction Impaired", "Road Contributing Circumstance",
"Ambient Light","Traffic Control Device","Improper Driving"),
column.labels=c("All Predictors","Backwards Elimination"),
add.lines=list(c("AIC", "11824.04", "11822.68"), c("BIC","11868.53","11854.45"), c("PRESS","4012.704","4011.383")),
type="text", out="stargazertable.txt", ci=TRUE, ci.level=0.05, single.row=TRUE)
## 
## Table 1. Regression Results
## ====================================================================================
##                                                 Dependent variable:                 
##                                -----------------------------------------------------
##                                                Injury Severity Level                
##                                      All Predictors         Backwards Elimination   
##                                           (1)                        (2)            
## ------------------------------------------------------------------------------------
## Pavement Friction Impaired      -0.042 (-0.045, -0.040)                             
## Road Contributing Circumstance  -0.081 (-0.085, -0.077)                             
## Ambient Light                  -0.171*** (-0.173, -0.169) -0.164*** (-0.166, -0.162)
## Traffic Control Device         -0.096*** (-0.098, -0.094) -0.097*** (-0.099, -0.096)
## Improper Driving                0.153*** (0.151, 0.155)    0.153*** (0.151, 0.155)  
## Constant                        2.598*** (2.596, 2.600)    2.582*** (2.580, 2.584)  
## ------------------------------------------------------------------------------------
## AIC                                     11824.04                   11822.68         
## BIC                                     11868.53                   11854.45         
## PRESS                                   4012.704                   4011.383         
## Observations                             4,253                      4,253           
## R2                                       0.014                      0.014           
## Adjusted R2                              0.013                      0.013           
## Residual Std. Error                0.971 (df = 4247)          0.971 (df = 4249)     
## F Statistic                     12.248*** (df = 5; 4247)   19.532*** (df = 3; 4249) 
## ====================================================================================
## Note:                                                    *p<0.1; **p<0.05; ***p<0.01
###Export stargazer as html too
stargazer(ols_levels, br_fit, title="Table 1. Regression Results",
dep.var.labels=c("Injury Severity Level"),
covariate.labels=c("Pavement Friction Impaired", "Road Contributing Circumstance",
"Ambient Light","Traffic Control Device","Improper Driving"),
column.labels=c("All Predictors","Backwards Elimination"),
add.lines=list(c("AIC", "11824.04", "11822.68"), c("BIC","11868.53","11854.45"), c("PRESS","4012.704","4011.383")),
type="html", out="stargazertable.html", ci=TRUE, ci.level=0.05, single.row=TRUE)
## 
## <table style="text-align:center"><caption><strong>Table 1. Regression Results</strong></caption>
## <tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="2"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="2" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td colspan="2">Injury Severity Level</td></tr>
## <tr><td style="text-align:left"></td><td>All Predictors</td><td>Backwards Elimination</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td></tr>
## <tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Pavement Friction Impaired</td><td>-0.042 (-0.045, -0.040)</td><td></td></tr>
## <tr><td style="text-align:left">Road Contributing Circumstance</td><td>-0.081 (-0.085, -0.077)</td><td></td></tr>
## <tr><td style="text-align:left">Ambient Light</td><td>-0.171<sup>***</sup> (-0.173, -0.169)</td><td>-0.164<sup>***</sup> (-0.166, -0.162)</td></tr>
## <tr><td style="text-align:left">Traffic Control Device</td><td>-0.096<sup>***</sup> (-0.098, -0.094)</td><td>-0.097<sup>***</sup> (-0.099, -0.096)</td></tr>
## <tr><td style="text-align:left">Improper Driving</td><td>0.153<sup>***</sup> (0.151, 0.155)</td><td>0.153<sup>***</sup> (0.151, 0.155)</td></tr>
## <tr><td style="text-align:left">Constant</td><td>2.598<sup>***</sup> (2.596, 2.600)</td><td>2.582<sup>***</sup> (2.580, 2.584)</td></tr>
## <tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">AIC</td><td>11824.04</td><td>11822.68</td></tr>
## <tr><td style="text-align:left">BIC</td><td>11868.53</td><td>11854.45</td></tr>
## <tr><td style="text-align:left">PRESS</td><td>4012.704</td><td>4011.383</td></tr>
## <tr><td style="text-align:left">Observations</td><td>4,253</td><td>4,253</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.014</td><td>0.014</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.013</td><td>0.013</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>0.971 (df = 4247)</td><td>0.971 (df = 4249)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>12.248<sup>***</sup> (df = 5; 4247)</td><td>19.532<sup>***</sup> (df = 3; 4249)</td></tr>
## <tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="2" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>